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Abstract. In a recent paper we proposed and compared various approaches to compute the ground 
state and dynamics of the Schrodinger-Poisson-Slater (SPS) system for general external potential 
and initial condition, concluding that the methods based on sine pseudospectral discretization in 
space are the best candidates. This note is concerned with the case that the external potential 
and initial condition are spherically symmetric. For the SPS system with spherical symmetry, 
via applying a proper change of variables into the reduced quasi-lD model we simplify the meth- 
ods proposed for the general 3D case such that both the memory and computational load are 
significantly reduced. 



1. Introduction 

The Schrodinger-Poisson-Slater (SPS) system, or the Schrodinger-Poisson-Xa system, serves as 
a local single particle approximation of the time-dependent Hartrcc-Fock equations. For its formal 
derivation we refer the readers to [T31 US] and references therein. The SPS system reads, in scaled 
form, 

(1.1) i9#(x,t) = --V 2 V' + V r e X t(x)V' + C F F F (x,^-a|V| i ^, xeR 3 , t> 0, 

(1.2) - V 2 V^(x,i) = |?M 2 , xel 3 , t>0 

(1.3) tp(x,t = 0) = tp (x), xel 3 . 

Here, the complex-valued ^(x, t) is the single particle wave function with limixi-^oo |^(x, t)\ = 
exponentially fast, 14xt( x ) is a given external potential, for example a confining potential, Vp(x, t) 
stands for the Hartree potential with decay condition limixi^oo Vp(x, t) = 0, and Cp (Cp > for 
repulsive interaction and Cp < for attractive interaction) and a (a > for electrons due to the 
physical nature) are interaction constants. 

There exit at least two important invariants of the SPS system (|1.1[) - (|1.3[) : the normalized mass 

(1.4) N^(;t)):=\\ip(;t)f = ( |^(x,t)| 2 dx-l, i > 0, 



and the energy 

(1.5) £(VM)):= 



\ |VV,| 2 + (Vext(x) + ^V P (x,t)\ H 2 - 



dx. t > 0. 



And, the ground state is defined as the minimizer of the following constraint minimization problem: 
Find 4> 9 ^S = {4>\ E{4>) < oo, N (4>) = 1} such that 

(1.6) E g :=E(<t> g )=minE(4>). 

<p£b 
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There is a series of analytical results on the SPS system in literatures; see, e.g., [H Q21 [TS] for 
its well-posedness and [SJ [7J [TTJ [TJ] for the existence and uniqueness of its ground states. Also, a 
detailed review was given in [TH] Section 1]. On the other hand, the numerics of the SPS system 
was considered in, e.g., [3j [6j El [16] . In particular, in [16] we proposed and compared different 
methods to compute the ground state and dynamics of the SPS system for general external potential 
and initial condition, ending with a conclusion that a backward Euler sine pseudospectral (BESP) 
method and a time-splitting sine pseudospectral (TSSP) method are the best choices to approximate 
the ground state and dynamics respectively. However, we have pointed out that when the external 
potential and initial conditions are with spherical symmetry, the original 3D problem reduces to a 
quasi-lD problem, for which the spectral-type methods BESP and TSSP cannot be directly extended 
and we suggested to apply the standard finite-difference to space derivatives (THl Remarks 3.2 and 
4.1]. The objective of this note is to propose spectral-type methods for the spherically symmetric 
case, which simplify the BESP and TSSP methods proposed in [16] for general 3D case, with the 
help of a proper change of variables for the reduced quasi-lD model. The simplified methods are still 
spectrally accurate in space, but reduce the memory cost from 0{ J 3 ) to 0{ J) and the computational 
cost per time step from 0(J 3 ln(J 3 )) to (9(Jln(J)), where J is the number of mesh nodes. 

The rest is organized as follows. In Section 2 we give a reduced quasi-lD model for the spherically 
symmetric case. In Section 3 simplified BESP and TSSP methods are proposed and in Section 4 
numerical results are reported. Finally, some concluding remarks are drawn in Section 5. 



2. A QUASI- ID MODEL REDUCED FROM SPHERICALLY SYMMETRIC SYSTEM 

Throughout this note, we assume that both the external potential V ex t and initial condition ipo 
are spherically symmetric, i.e., K> x t(x) = V cx t(f) and V'o(x) = ipo(r) with r = |x|. In this case, the 
solution -0 of (|l.ip - (|1.3[) and the ground state <j> g are also spherically symmetric, i.e., 

V>(x, t) = ip(r, t), 9 (x) = 4>{r), x e M 3 , t>0. 

Thus, the SPS system (|1.1[) - (|1.3I) collapses the following quasi-lD problem 

1 d ( 2 dtp 



(2.1) id t ^(r,t) = - — —^r 2 -^j+V cxt (r)ijj + CpVp(r,t)i;-a\ib\ 3 TP, < r < oo, t > 0, 

<"> -jMfO* 2 ^) -1 * 1 '' 0<r< ~' '- - 

(2.3) VO, t = 0) = Vo(r)> < r < cx), 
with boundary conditions 

(2.4) d r ip(0,t) = d r Vp(0,t)=0, lira. ip(r,t) = 0, lim rV P (r, t) = — , t > 0, 

due to the decay conditions of ip and Vp , and the Green function of the Laplacian on R 3 [10] . 
Introducing 

(2.5) U{r,t) = 2y/Trrip(r,t), V(r, t) = 4irrVp(r, t), < r < oo, t > 0, 
a simple computation shows 

(2.6) i.* (r^) = 1 3 rrU 11 = ^ drrV . 
r 2 dr \ dr J rr ' r 2 dr \ dr J 4irr rr 
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We remark here that the similar technique has been used in [TJIH]. Plugging the above into (|2.1j) - (l2.4l) . 
we obtain 

(2.7) id t U{r,t) = -\d rr U + V cxt {r)U + ^-V{r,t)U-a{2^rY^\U\^U, < r < oo, t > 0, 

2 inr v /ii 

(2.8) - d rr V(r, t) = - \U\ 2 , < r < oo, t > 0, 

r 

(2.9) W(r,i = 0) =W (r) =2V^rrV>o(r), < r < oo, 

(2.10) W(0, t) = V(0, f) = 0, lim W(r, t) = 0, lim V(r, t) = 1, t > 0. 

Also, the above problem conserves the mass 

/•OO 

Af(U{;t)) := ||W(-,i)|| 2 = / |WM)| 2 dr = JV(^(-,t)) = l, t > 0, 



and the energy 



f («(•,*)) := 



~ \d r U\ 2 + (w»0 + ^V(r,t)) |W| 2 - ^ (2^) 1 



dr = E(ip(-,t)), t>0. 



In what follows we will take the problem (|2.7[) - (|2.10[) as the starting model to propose efficient 
numerical methods. After we get the solution U of (12.7[) - (|2.10[) . the solution -0 of (|2.1|l - (|2.4|l is 
obtained as 

(2 11) Mr t) = — ( W ( r '*)/ r ' r>0 ' t >0 

{ ] n , ) 2^{ d r U(r,t) =\im s ^ 0+ U(s,t)/s, r = 0, - u ' 



Meanwhile, the minimization problem (|1.6|) to define the ground state collapses to: 
Find (fig e S = {up | £(<p) < oo, Af(<p) = 1, <^>(0) = 0} such that 

(2.12) £„ := £(<p g ) = mm£(ip). 

Again, after we get the minimizer of ()2. 12|) . the ground state <j> g of (|2.1|) - (j2.3[) is obtained as 

1 f l Pg( r )/ r , r >0, 

2y/TT \ d r ip g {r) = lim s _> + <p g (s)/s, r = 0. 

3. Efficient numerical methods 



(2.13) <f> g (r) 



3.1. Backward Euler sine pseudospectral method for ground state. Choose a time step 
At > and set t n = nAt for n = 0, 1, . . . . Similar as |161 Section 2], for the minimization problem 
(|2.12[) . we construct the following gradient flow with discrete normalization (GFDN): 

(3.1) 

d t ip(r,t) = \d rr ip - V CKt (r)p - -r^-V(r,t)<p + a(2y / nr) 3 \tp\ 3 <p, < r < oo, t n < t < t n+ i, 
(3.2) 

- d rr V(r,t) = - \tp\ 2 , < r < oo, t > 0, tp(r,t n+1 ) := <p(r,t+ +1 ) = ' , n > 0, 

r ||W*„+i)|| 

(3.3) 

y>(r,i = 0) = ^o(r), < r < oo, with N(tp Q ) = l, 
(3.4) 

<^(0, i) = V(0, t) = 0, lim <p(r, t) = 0, lim V(r, t) = 1, i > 0, 

r— >oo r— >oo 

where ^?(r, t^) '■— Iina t _^ t ± tp(r,t) for < r < oo. In practical computation, we truncate the above 
problem into an interval [0, R] with R > sufficiently large, together with Dirichlet boundary 
conditions 

<p(0, t) = <p(R, t) = V(0, t) = 0, V(i?, t) = 1, t>0. 
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Introducing a linear translation (homogenization) V(r, t) = V(r, t) — r/R for < r < R, 



(3.5) - d rr V(r,t) = -d rr V(r,t) = - Id 2 , < r < R. V(0, t) = V(R, t) = 0, t>0. 

r 

Then we discretize the problem in space by sine pseudospectral method and in time by a backward 
Euler integration similar as that used in [16] . Choose a mesh size h r = Ar — R/J with some 
even integer J > 0, and denote the grid points as Tj = jh r for j = 0, 1, . . . , J. Let tp" ps tp(rj,t n ) 

and V™ fa V(rj,t„), and denote p" = |<^™| /Vj. Choosing </?9 = (po(fj), a backward Euler sine 
pseudospectral discretization (BESP) reads: for n — 0, 1, . . . , 



(3.6) 

<Pj-<P? 
At 

(3.7) 
(3.8) 



(^ + )L- 



Wri) + ^ + ^-a(W 



j = l,2,...,J-l, 



P", J = 1,2,...,J-1, ^+ = ^+ = Vo=V3 



n+l _ 



, j = 1,2,..., J- 1, with 



j-i 

3=1 



where _D£ r is the sine pseudospectral approximation of <9 rr , defined via 



J-i 



(3-9) - (£>")|. = 5> 2 & (^) fe sin 

fc=i \ 17 

with the discrete sine transform coefficients 



, j = 1,2,..., J- 1, 



J-i 



(3.10) 



3 = 1 



j'fcvT 



■J 



, , z v — ,„ / yre/i \ re/i 



R 



Similar as [IB], the linear system (|3.6l) - (|3.8p can be iteratively solved efficiently in phase space with 
the help of discrete sine transform. After we get the stationary solution (<p g )j of the above problem, 
the ground state (</> g )j ~ <t>g{ r j) °f (I2.1[) - (|2.3[) is achieved via 



(3.11) 



1 



,J gio 



J = 1,2,..., J, 



2 V* I Efe=l t*k(<Pg) k , j = 0. 



Note that the above numerical method is spectrally accurate and it works only when V cx t is spherically 
symmetric. Compared with the pseudospectral method proposed in |16| for general 3D problem, the 
memory cost is reduced from 0(J 3 ) to O(J) and computational cost per time step is reduced from 
0(J 3 ln(J 3 )) to 0(Jln(J)). 



3.2. Time-splitting sine pseudospectral method for dynamics. Again, we truncate the prob- 
lem into an interval [0, R], and introduce the linear translation Q3.5P for V into (|2.7|) - (|2.10|) such that 
both U and V satisfy homogeneous Dirichlet boundary conditions. Similar as [16] . for computing the 
dynamics, we first apply the time-splitting technique to decouple the nonlinearity and then use sine 
pseudospectral method to discretize the spatial derivatives. Let ps U(rj,t n ) and Vj ~ V{rj,t n ). 
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Choose 1$ = Uo(rj), a second-order time-splitting sine pseudospectral (TSSP) discretization reads: 



j-i 



(3.12) 


Uf> = 


H ex p{- 






k— 1 


(3.13) 


Uf = 


exp /— iAt 






.7-1 


(3.14) 




= £exp{- 






fe=i 


for n > 0, and j = 


1,2,..., J- 



jklt 

Cp ^(i) , c t 



-V 



— a 



U 



(i) 



jkn\ 

j r 



Here, V,-^ is obtained 
pseudospectral method (similar as [16j Section 3.2]), i.e. 

J-i . _ 



(3.15) vf^EMrV 1 ^ 



k=l 



jkn 



with p 



(1) 



U 



(i) 



j = l,2....,J-l. 



Again, after we get the solution W from (|3~12t - (|3~l"5)) . the solution ipf w ip(rj,t n ) of ((2~T1) - ((2~3 



is achieved via 
(3.16) 



j = l,2,...,J, 



1 ELiw(w n ) fc , i = o. 



The above method is explicit, spectrally accurate in space and second-order accurate in time and 
it works only when both T4 x t and ipo are spherically symmetric. Again, compared with the method 
proposed in [16] for general 3D problem, the memory cost is reduced from 0(J 3 ) to O(J) and 
computational cost per time step is reduced from 0(J 3 ln(J 3 )) to 0(Jln(J)). In addition, similar 
as [TB], we have, 

Lemma 3.1. The TSSP method i3.1ty) - [3.15\) is normalization conservation, i.e., 

J-i J-i 



so it is unconditionally stable in L 2 -norm. 



4. Numerical results 

Numerical results are reported in this section to demonstrate the accuracy and efficiency of the 
proposed methods, and we choose V ex t = \r 2 '', Cp = 100 and a = 1 in (|2.ip as the example. 
For computing the ground state, the "exact" solution 4> g (benchmark) is achieved by applying a 
backward Euler finite-difference (BEFD) discretization to a GFDN of the quasi-lD model (12.1[) - (|2.3p 
with Dirichlet boundary conditions of 4> and Robin boundary conditions of Vp [16j . 4> g is computed 
in a ball < r < 8 with a very fine mesh size h r = 1/64. Let (j>g be the approximations obtained 
from BESP method ([3?6])-(j3l]), Fig. [T] plots <j> g and <j>* with h r = 1/2, and the error \(j> g - $j\ with 
different h r . The results show that the BESP method (|3.6p - (|3.8|) gives the approximation of ground 
state with spectral order of accuracy in space; and therefore, it is more efficient in implementation 
than the standard finite-difference discretization for spherically symmetric case and the spectral- 
type method proposed in [16] for general 3D case. Similar accuracy and efficiency conclusions can 
be drawn for TSSP method (|3~T2> (|3~T4l) . Fig. H plots the evolution of \ip n \ for < t n < 10 when 
ipo = (27r) 3 / 4 exp (— r 2 /4). Here, the computation is carried out in a ball < r < 16, with h r = 1/16 
and At = 0.01. 
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Figure 1. Accuracy analysis for BESP method: (1) (f> g obtained from BEFD 
method with h r = 1/64 as benchmark and 0^ obtained from BESP method with 



h r = 1/2 (left figure); (2) error 



bg \ with different h r (right figure) 



0.05 




Figure 2. Dynamics computed by TSSP method: evolution of |?/? n | up to time t n = 10. 



5. Concluding remarks 

In this study we considered the numerics of the spherically symmetric SPS system and simplified 
the spectral-type methods proposed in our recent paper [TB, to compute the ground state and dynam- 
ics for general external potential and initial condition. The simplification is achieved by introducing 
a proper change of variables into the reduced quasi-lD model. The simplified methods still admit 
spectral order of accuracy in space, with significantly less demand on memory and computational 
load, and is more efficient in implementation than the standard finite-difference approaches for the 
spherically symmetric case. Note that the simplified methods only work for the system with spherical 
symmetry, and for the general case we still suggest to apply the methods proposed in [16j . Also, 
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the results in this study are applicable to the Schrodinger-Poisson and Schrodinger-Newton systems 
(a = in ([TT ]l -([r5 ]) ) as well. 
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